*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0

* This scripts defines various programmes in Stata programming syntax that will
* be called by other Stata scripts. 

* GRID ************************************************************************* 
* Program creating a 70 x 70 km grid that consists of 250 x 250 m cells

	capture program drop GRID
	program define GRID
	clear
	set obs 280 // 70 km grid length / 250 m cell size // 400 gives 100 km grid length
	gen x=_n 
	gen y = x
	fillin y x 
	gen double grid_x = x*250
	gen double grid_y = y*250
	sum grid_x
	replace grid_x = grid_x - r(mean) + round(`1',250)
	sum grid_y
	replace grid_y = grid_y - r(mean) + round(`2',250)
	drop x  y _fillin
	end

* AGGCLUSTER2 ******************************************************************
* Program that operationalizes lines 18-29 of Aglorithm 1
* Takes clusters generated in C++ as inputs
* Returns aggregated clusters of sufficient size as outputs (prime locations) 

	capture program drop AGGCLUSTER2
	program define AGGCLUSTER2 	// Syntax uses the following arguments
								// 1 metro_id 
								// 2 distance threshold in meters
	preserve
	qui keep if metro_id == `1'		// keep metro
	display "...working on aggregatig clusters in metro `1'..."
	qui collapse (first) TrimmedClusterX TrimmedClusterY Cemp, by(TrimmedClusterID metro_id)	// collapse data set to trimmed clusters
	qui drop if TrimmedClusterID == .	// only keep positive observations
	qui drop if TrimmedClusterID == 0
	qui egen TrimmedClusterNID = group(TrimmedClusterID) // generate temporary continuous ID
	qui gen grid_y = TrimmedClusterY					// Generate mercator y coordinate for projection correction
	qui gen NewTrimmedClusterID = .		// new variable for new allocation to PL IDs
	* Gen distance matrix
	qui sum TrimmedClusterNID
	qui local Ncluster =r(max)
	qui foreach clusterNID of numlist 1/`Ncluster' {
		qui sum TrimmedClusterID if TrimmedClusterNID == `clusterNID', meanonly // get cluster ID 
		local cluster = int(r(mean))
		qui sum TrimmedClusterX if TrimmedClusterID == `cluster' // get X-coordinate of cluster
		local X = (r(mean))
		qui sum TrimmedClusterY if TrimmedClusterID == `cluster' // get Y-coordinate of cluster
		local Y = (r(mean))
		qui sum Cemp if TrimmedClusterID == `cluster' // get emp of cluster
		local E = (r(mean))
		qui gen  dist`cluster' = $bmd *sqrt([TrimmedClusterX-[`X']]^2+[TrimmedClusterY-[`Y']]^2) + int(runiform()*1000)/1000 // We add a microscopic distance value to avoid later ties when evaluation min(dist) in a couple of lines
		qui gen Cemp_pair`cluster' = [`E']
	}
	* reshape to long
	qui reshape long dist Cemp_pair , i(TrimmedClusterID) j(TrimmedClusterID_pair)
	qui egen temp = min(dist) if TrimmedClusterID != TrimmedClusterID_pair, by(TrimmedClusterID)		  // identify potential pair
	qui gen pair = temp == dist							  // identify potential pair	
	qui replace pair = 0 if temp > `2'					  // identify potential pair
	qui replace NewTrimmedClusterID = TrimmedClusterID_pair if pair == 1 &  Cemp_pair > Cemp // update ID if pair has larger employment
 	qui collapse  NewTrimmedClusterID , by(TrimmedClusterID metro_id)  			// gen mapping from old to new ID
	qui replace NewTrimmedClusterID = TrimmedClusterID if NewTrimmedClusterID == .  // gen mapping from old to new ID
	 qui keep TrimmedClusterID NewTrimmedClusterID metro_id // save key variables
	capture mkdir "$temp/PLtemp"								// Create temporary directory	 
 	qui save "$temp/PLtemp/NewClusterIDtemp.dta", replace	
	restore 
	 qui merge m:1 metro_id TrimmedClusterID using "$temp/PLtemp/NewClusterIDtemp.dta" // merge output to baseline data
	 qui drop _m
	qui replace TrimmedClusterID = NewTrimmedClusterID if metro_id == `1'		// replace cluster ID with new ID
	drop NewTrimmedClusterID	
	// Update coordinates
	qui egen tempX = mean(grid_x) if TrimmedCluster > 0 & metro_id == `1', by(TrimmedClusterID metro_id)
	qui egen tempY = mean(grid_y) if TrimmedCluster > 0 & metro_id == `1', by(TrimmedClusterID metro_id)
	qui egen tempCemp = mean(emp) if TrimmedCluster > 0 & metro_id == `1', by(TrimmedClusterID metro_id)
	qui replace TrimmedClusterX = tempX if metro_id == `1'
	qui replace TrimmedClusterY = tempY if metro_id == `1'	
	qui replace Cemp = tempCemp if metro_id == `1'
	qui drop tempX tempY tempCemp
	preserve
		qui filelist, dir("$temp/PLtemp")
		qui levelsof filename, local(files)
		foreach f of local files{
		qui 	erase "$temp/PLtemp/`f'"
		}
	restore 	
end 	

* CM ***************************************************************************
* Program to generate a minimum cluster mass around microscopic clusters
capture program drop CM
program define CM 		// Syntax uses the following arguments
						// 1 metro_id 
						// 2 cluster ID 
						// 3 fill X around median coordinate
	display "...buildig mass for cluster `2' in metro `1'..."
	if `3' > 0 {
		qui sum grid_x if  metro_id == `1' & clusterID == `2', d
			local medX = r(p50)
		qui sum grid_y if metro_id == `1' & clusterID == `2', d
			local medY = r(p50)
		qui replace clusterID = `2' if grid_y>=`medY'-`3'*250 &  grid_y<=`medY'+`3'*250 & grid_x>=`medX'-`3'*250 &  grid_x<=`medX'+`3'*250
		
	}
	else {
	}
end
	
* CH ***************************************************************************
* Program that operationalizes line 30 of Algorithm 1
* Takes aggregated clusters generated by AGGCLUSTER2 as input
* Creates a convext hull by filling gaps in horizontal or vertical dimension	
	
capture program drop CH
program define CH 		// Syntax uses the following arguments: 
						// 1 metro_id 
						// 2 PL ID 
						// 3 Order of rook contiguity cells around median coordinate that are considered part of the prime location
	display "...buildig convex hull for PL `2' in metro `1'..."
	if `3' > 0 {
		display "...adding center to PL `2'"
		qui sum grid_x if  metro_id == `1' & PLID == `2', d
			local medX = r(p50)
		qui sum grid_y if metro_id == `1' & PLID == `2', d
			local medY = r(p50)
		qui replace PLID = `2' if grid_y>`medY'-`3'*250 &  grid_y<`medY'+`3'*250 & grid_x>`medX'-`3'*250 &  grid_x<`medX'+`3'*250
		
	}
	else {
	}
	
	* fill Xs
	qui sum grid_y if metro_id == `1' & PLID == `2'
		local ymin = r(min)
		local ymax = r(max)
	qui local i = `ymin'
	qui while `i' <= `ymax' {
		sum grid_x if  metro_id == `1' & PLID == `2' & grid_y == `i'
		local xmin = r(min)
		local xmax = r(max)
		replace PLID = `2' if  metro_id == `1' & grid_y == `i' &  grid_x >= `xmin'  & grid_x <= `xmax'
		local i = `i'+250
	}
	* fill Ys
	qui sum grid_x if metro_id == `1' & PLID == `2'
		local xmin = r(min)
		local xmax = r(max)
	qui local i = `xmin'
	qui while `i' <= `xmax' {
		sum grid_y if  metro_id == `1' & PLID == `2' & grid_x == `i'
		local ymin = r(min)
		local ymax = r(max)
		replace PLID = `2' if  metro_id == `1' & grid_x == `i' &  grid_y >= `ymin'  & grid_y <= `ymax'
		local i = `i'+250
	}	
end	

* SPLITTER ***************************************************************************
* Implements the border-distance condition in line 25 of Algorithm 1
* If the between the borders of to aggregated clusters is too large, they will be treated as separate PLs
* Takes aggregated clusters generated by AGGCLUSTER2 as input
capture program drop SPLITTER 
program define SPLITTER // Syntax: 1 cutoff distance
	local 1 = `1'/1000 	// Rescale split distance from m to km
 	qui capture mkdir "$temp/TEMPSPLIT" // Generate temporary working directory
	preserve
		qui filelist, dir("$temp/TEMPSPLIT")
		qui levelsof filename, local(files)
		qui foreach f of local files{
		qui	erase "$temp/TEMPSPLIT/`f'"
		}
	restore 
	qui save  "$temp/TEMPSPLIT/plworking", replace // save data set
	qui drop if PLID==. // drop non-PL cells
	qui levelsof metro_id	, local(MID)
	qui sum PLID
	local max_PLID=r(max)
	foreach m_ID of local MID { 												// begin outer loop over metros
		display "...considering PLs in metro `m_ID' for splitting..."
		u "$temp/TEMPSPLIT/plworking", clear										// start from full data set
		qui drop if PLID==. // drop non-PL cells
	 	qui keep if metro_id==`m_ID'	
		* count PLS within metro_id 
		levelsof PLID		, local(N_PLS)
		* di `N_PLS'
	foreach PL_i of local N_PLS {												// begin inner loop over PLs wihtin metro
		qui u "$temp/TEMPSPLIT/plworking", clear
		qui drop if PLID==. // drop non-PL cells
		qui keep if metro_id==`m_ID'
		qui keep if PLID==`PL_i' 
		display "...PLs `PL_i'..."
	* identify gap in PL 
	* overlay a square with min and max of lat lon
						qui format grid_x  %12.0g 
						capture format grid_y long
						*replace grid_x=-1
						qui foreach r in x y{
								sum grid_`r'
								local `r'_min=r(min)
								local `r'_max=r(max)
							}	
						* identify missing X-cells 
						qui forvalues r=`x_min'(250)`x_max' {
								*di `r'
								local  k= abs(`r')			    
								gen  cell_x_`k'=`r'  if `r'==grid_x
								egen cell_x_`k'_tot=sum(cell_x_`k')
								drop cell_x_`k'
								sum cell_x_`k'_tot 
								if r(mean)!=0 drop cell_x_`k'_tot
								capture rename cell_x_`k'_tot E_x_`k'
								capture replace  E_x_`k'=`r'
								capture format  E_x_`k' %12.0g
								}
						* identify missing Y-cells 
						qui forvalues r=`y_min'(250)`y_max' {
								*di `r'
								local  k= abs(`r')			    
								gen  cell_y_`k'=`r'  if `r'==grid_y
								egen cell_y_`k'_tot=sum(cell_y_`k')
								drop cell_y_`k'
								sum cell_y_`k'_tot 
								if r(mean)!=0 drop cell_y_`k'_tot
								capture rename cell_y_`k'_tot E_y_`k'
								capture replace  E_y_`k'=`r'
								capture format  E_y_`k' %12.0g
								}
						* group the two sub-cluster in a given PL 
						capture egen median_Empty_y=rowmedian(E_y_*)
						capture egen median_Empty_x=rowmedian(E_x_*)
						qui gen group_x=0
						qui gen group_y=0
						qui foreach r in x y {
							capture		replace group_`r'=1 if grid_`r'<median_Empty_`r'
							}
						qui egen xygroup=group(group_x group_y)
			* Now calculate shortest edge-to-edge distance
			* export file for geonear for the available groups 
			qui levelsof xygroup, local(groups)
			qui  sum xygroup 
			qui if r(sd)!=0{
			* insert if condition to only calculate if there non-adjacent PL-islands	
						qui foreach x of local groups {
							preserve 
								qui keep if xygroup==`x'
								qui keep square_id lat lon 
								rename square_id nborid
								rename lat nborlat
								rename lon nborlon
								save  "$temp/TEMPSPLIT/nbordata_`x'", replace 
							restore 
						}
						qui foreach x of local groups { 
							qui geonear  square_id lat lon using "$temp/TEMPSPLIT/nbordata_`x'" , genstub(g`x'_) n(nborid nborlat nborlon)
							qui replace km_to_g`x'_=. if   km_to_g`x'==0
							qui egen min_dist_`x'=min(km_to_g`x')
						}
					* break up if parameter for min distance is surpassed and rename PL_ID 
					qui gen PLID_sub="a"  if min_dist_1>=`1' & xygroup==1
					qui replace PLID_sub="b"  if min_dist_1>=`1' & xygroup==2
					qui egen PLID_new=concat (PLID PLID_sub), punct(_)
			}
		qui sum xygroup 
		qui if r(sd)==0{
					tostring PLID, g(PLID_new)
		}					
		qui drop if PLID==. /// cleans up non-PL cells 
		
		*local PLID_new=PLID_new[1]
 		qui save "$temp/TEMPSPLIT/IC_`m_ID'_`PL_i'" ,replace 
 	}																			// close inner PL loop
	}																			// close outer metro loop
	* assign new IDs 
	display "...process new PL IDs..."
			qui clear 
			qui foreach m_ID of local MID {
				clear 
				qui forvalues r=1/`max_PLID'{
				 capture append using "$temp/TEMPSPLIT/IC_`m_ID'_`r'"
				}
				  
				qui drop PLID  
				qui egen PLID=group(PLID_new)
				qui drop PLID_new
				qui save "$temp/TEMPSPLIT/PLS_`m_ID'", replace 
		}
	* merge all data 	
			qui clear 
			qui foreach m_ID of local MID {
				 append using "$temp/TEMPSPLIT/PLS_`m_ID'"
				 }
			qui save "$temp/TEMPSPLIT/splitPLs.dta", replace
	display "...update PL IDs..."
	qui u "$temp/TEMPSPLIT/plworking" , clear
	qui drop PLID
	qui merge 1:1 grid_x grid_y  using 	"$temp/TEMPSPLIT/splitPLs.dta", keepusing(PLID)	
	qui 	drop _m	
	preserve
		qui filelist, dir("$temp/TEMPSPLIT")
		qui levelsof filename, local(files)
		qui foreach f of local files{
		qui	erase "$temp/TEMPSPLIT/`f'"
		}
	restore 
end 


* DISCONNECT *******************************************************************
* Program that operationalizes lines 30-34 of Alogorithm 1
* Disconnects prime locations that are very large and have a strange geography (thin waist)
* Takes aggregated clusters generated by AGGCLUSTER2 as input

capture program drop DISCONNECT
program define DISCONNECT // Syntax uses the following arguments: 
						  // 1: prime location id
					      // 2: Size threshold
	* Gen a flag variable to keep track of disconnections
		display "...evaluating size of prime location `1'..."
		capture gen flag_DCPL = 0
		* Disconnect horizontally
				qui sum PLID
					local PLIDmax = r(max)										// Save the number of prime locations
					qui sum grid_y  if PLID == `1'								// Measure vertical size
					local y_size = r(max)-r(min)
					local y_min = r(min)
					local y_max = r(max)
				if `y_size' < `2' {												// PL is small enough
					* do nothing
					}
				else {
					* Check where thinnest
					display "...evaluating potential disconnection..."
					qui egen x_min = min(grid_x) if PL == `1', by(grid_y)			// Get horizontal size by row
					qui egen x_max = max(grid_x) if PL == `1', by(grid_y)
					qui gen x_size =  x_max - x_min if x_max - x_min <= 250
					qui sum grid_y  if PLID == `1' , d								// Find vertical center of PL
					qui gen y_centdist = abs(grid_y - r(p50)) if x_size != .		// Gen vertical distance from center	
					qui egen y_centdist_min = min(y_centdist) if x_size != .		// Find smallest vertical distance from center
					* Check if criteria for disconnecting are met
						qui gen CHECK = x_size <= 250 & y_centdist == y_centdist_min 	/// Must be thick, closest to center, and close to center
						& y_centdist <= `2'/2 ///
						& grid_y <= (`y_max'-1000) & grid_y >= (`y_min'+1000)		// Do not split to close to edges
						qui sum CHECK
						local CHECKmax = r(max)
					if `CHECKmax' == 0 {
						display "...prime location NOT disconnected..."
					}
					else {
						display "...prime location is being disconnected..."
						qui sum grid_y if x_size <= 250 & y_centdist == y_centdist_min ///
						& y_centdist <= `2'/2, meanonly								// Find y-value of small row that is closest to the center conditional on being close
						qui replace PLID = . if grid_y == r(mean)					// If central, cut at closest narrow row
						qui replace PLID = `PLIDmax'+1 if grid_y <  r(mean) & ///
						PLID == `1'													// Update PLID
						qui replace flag_DCPL = 1 if PLID == `PLIDmax'+1			// Flag PL as disconnected
					}
					qui drop CHECK													// Drop marker for diconnection criterion
				}
				* Disconnect vertically
				qui sum PLID
					local PLIDmax = r(max)											// Save the number of prime locations
					qui sum grid_x  if PLID == `1'									// Measure horizontal size
					local x_size = r(max)-r(min)	
					local x_min = r(min)
					local x_max = r(max)
					if `x_size' < `2' {												// PL is small enough
					* do nothing
					}
				else {
					* cut where thinnest
					display "...evaluating potential disconnection..."
					qui egen y_min = min(grid_y) if PL == `1', by(grid_x)			// Get vertical size by column
					qui egen y_max = max(grid_y) if PL == `1', by(grid_x)
					qui gen y_size =  y_max - y_min if y_max - y_min <= 250
					qui sum grid_x  if PLID == `1' , d								// Find horizontal center of PL
					qui gen x_centdist = abs(grid_x - r(p50)) if y_size != .		// Gen horizontal distance from center	
					qui egen x_centdist_min = min(x_centdist) if y_size != .		// Find smallest horizontal distance from center
					* Check if criteria for disconnecting are met
						qui gen CHECK = y_size <= 250 & x_centdist == x_centdist_min /// Must be thick, closest to center, and close to center
						& x_centdist <= `2'/2										///
						& grid_x <= (`x_max'-1000) & grid_x >= (`x_min'+1000)		// Do not split to close to edges

						qui sum CHECK
						local CHECKmax = r(max)					
					if `CHECKmax' == 0 {
						display "...prime location NOT disconnected..."
					}	
					else {
						display "...prime location is being disconnected..."
					qui sum grid_x if  y_size <= 250 & x_centdist == x_centdist_min ///
					& x_centdist <= `2'/2, meanonly									// Find y-value of small row that is closest to the center conditional on being close
					qui replace PLID = . if grid_x == r(mean)						// If central, cut at closest narrow row
					qui replace PLID = `PLIDmax'+1 if grid_x <  r(mean) & PLID == `1'	// Update PLID
					qui replace flag_DCPL = 1 if PLID == `PLIDmax'+1				// Flag PL as disconnected
					}
					qui drop CHECK													// Drop marker for diconnection criterion
				}				
				qui foreach tempvar in x_min x_max x_size y_centdist y_centdist_min y_min y_max y_size x_centdist x_centdist_min {
					capture drop `tempvar'
				}
			end	

* BDD **************************************************************************
* Program that generates the boundary distance trends displayed in Figure B.1.1

capture program drop BDD
program define BDD	// Syntax
					// 1 List of outcome variables in ""
					// 2 Distance from PL boundary variable  
					// 3 p-value in ""
					// 4 Bin size in km
					// 5 Outer limite of gradient (add 1km extra as one bin has to be ommitted)
					// 6 Inner limit of gradient
	* Preserve data set
		preserve
	* Read user choices of outcomes
		local valvarlist = "`1'" 
	* Read user choices of p-value
		local pvalue = "`3'"		
	* Read user choice of PL distance variable
		qui gen distance_toPLborder = `2'
	* Read user choice of bin size
		local binsize_1km= `4'
	* Read user choice of outer limit
		local max_x_axis= `5'
	*  Read user choice of inner limit 
		local drop_inner=-`6'
	 	 
	display "..working on p-value `3'..."	 
		 
	* generate bins (flooring and ceiling sorts them into distinct bin (no overlapping 0 bin))
		qui g DISTB=floor(distance_toPLborder*10)/10 if distance_toPLborder<0
		qui replace DISTB=ceil(distance_toPLborder*10)/10 if distance_toPLborder>=0
		qui forvalues r=-2(`binsize_1km')-.2{
		* di `r'
		replace DISTB=`r'  if distance_toPLborder>=`r' &distance_toPLborder<`r'+`binsize_1km' &	distance_toPLborder<0
		}
		qui forvalues r=.2(`binsize_1km')1{
		* di `r'
		replace DISTB=`r'  if distance_toPLborder>=`r'-`binsize_1km'  &distance_toPLborder<`r' &	distance_toPLborder<1
		}
		qui replace DISTB=floor(DISTB) if distance_toPLborder>1
	 * dropping observations deep within PL for lack of sufficient observations (add 10 meters to avoid precision problems)
		qui drop if distance_toPLborder<`drop_inner'
		qui drop if distance_toPLborder>=`max_x_axis'
	
	* Gen dummies for bins
		qui tab DISTB  , gen(PLDB)
		local N_bins=r(r)
		local N_bins_less_ommited=`N_bins'-1
	* Bins inside PPL
		qui tab DISTB  if distance_toPLborder<1
		local binsunder1km=r(r)
		local binsover1km=`binsunder1km'+1
	* Label bins
		qui forval bin = 1/`binsunder1km' {
		local lower =   (`bin'-1)*`binsize_1km'+`drop_inner'
		local upper = (`bin')*`binsize_1km'+`drop_inner'  
		label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"  
		}		
		qui forval bin = `binsover1km'/`N_bins' {
		local lower =   (`bin'-`binsunder1km')
		local upper = 	(`bin'-`binsunder1km')+1  
		label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"
		}
	* RUN BIN REGRESSIONS WITH METRO FIXED EFFECT
		qui foreach var of varlist `valvarlist' {
		eststo: areg `var' PLDB1-PLDB`N_bins_less_ommited' if developable == 1  , abs(metro_id) cluster(metro_id)
				estadd scalar metro_areas = e(N_clust)
				estadd local metro_FE = "Yes"
				forval bin = 1/`N_bins_less_ommited' {
					scalar b`var'_`bin' = _b[PLDB`bin'] 
					scalar se`var'_`bin' = _se[PLDB`bin'] 
					}
				}
			qui esttab using "$temp/TAB-PLgradients_3USCITIES'.rtf", replace b(2) se(2) onecell label compress stats(metro_areas metro_FE N r2, fmt(%18.3g )) drop(_cons ) ///
			title ("Distance from  prime location gradients: Distance grid cell models") modelwidth(6) nogap star(* 0.1 ** 0.05 *** 0.01)  ///
			addnote("Unit of observation are 250 x 250 m grid cells (excluding undeveloped cells" "Standard errors clustered on metro areas.")
			eststo clear			
	* Read results of regressions for plotting
		clear 
		qui set obs `N_bins_less_ommited'
		qui gen bin = _n
		qui foreach var in `valvarlist' {
			gen b`var' 	= .
			gen se`var' = .
			}
	* READ SCALARS 
		qui foreach var in `valvarlist' {
		qui forval bin = 1/`N_bins_less_ommited' {
			replace  b`var' =  b`var'_`bin' 	if bin == `bin' 
			replace  se`var' =  se`var'_`bin' 	if bin == `bin'
			}
		}
	* COMPUTE CIs
		qui foreach var in `valvarlist' {
		gen cup`var'  = b`var'  + 1.645* se`var' 
		gen clo`var'  = b`var'  - 1.645* se`var' 
		}
					
 	* GEN DIST
		qui gen dist= ((bin-1)*`binsize_1km'+`binsize_1km'/2)-(-`drop_inner') if bin<=`binsunder1km'
		qui replace dist=bin-`binsunder1km'+.5 if bin>`binsunder1km'
		qui gen  ln_dist=log(1+dist)
	
	* Loop pver outcomes
			qui foreach var in `valvarlist' {
				* Read p value for figure title
				local title= subinstr("`3'","_",".",1)
	
			* Recover outer border effect for labelling
				qui sum ln_dist if dist>0
				local mindist = r(min)
				qui sum b`var' if dist > 0
				local OBE = round(r(max), 0.01)
				local OBE_ft : display %3.2f `OBE'
			* Plot 
			twoway /// (rarea  cup`var'_`prime' clo`var'_`prime' dist, sort color(gs14) ) (line b`var'_`prime' dist, lcolor(black))   ///
			 (rarea  cup`var' clo`var' ln_dist if dist<0, sort color(red%25) lcolor(red%0)) (rarea  cup`var' clo`var' ln_dist if dist>0, sort color(red%25) lcolor(red%0)) ///
			 (line b`var' ln_dist if dist<0, lpattern(shortdash) lcolor(red)) (line b`var' ln_dist  if dist>0, lpattern(shortdash) lcolor(red)) 	///
					,  xlabel(-2.995 "0.95" -2.302 "-0.9" -1.38629 "-0.75" -0.693 "-.5"  0 "0" 0.693 "1" 1.253 "2.5" 1.79 "5" 2.34 "10"  3.045 "20") ylabel(#5) graphregion(color(white) lwidth(thick))  yline(0) xline(0)  /// xscale(log) xlabel(-2.302 "-0.9" -1.38629 "-0.75" -0.693 "-.5" -0.287 "-.25" 0 "0" 0.693 "1" 1.253 "2.5" 1.79 "5" 2.34 "10"  3.045 "20")
			xtitle("Distance from prime location (km)")  title("Cutoff at `title' %",  size(medium)) name(`var', replace) legend(off) note("Outer border effect: `OBE_ft'")
			qui graph save "$temp/_FIG_`var'_p`pvalue'.gph", replace		
				
			} // End loop over outcomes
		// Restore variables
		restore
	// End loop 
	end
	
* Script ends